#include <iostream>
#include "RestrictionOperatorFactory.h"
#include "InterpolationOperatorFactory.h"
#include "MultigridSolverFactory.h"
#include "VectorFunctionFactory.h"
#include <iomanip>
#include <fstream>
#include <sys/types.h>
#include <sys/stat.h>
#include <lapacke.h>

int main()
{
	char c;
	size_t DIM;
	std::string MultigridSolver_Id;
	std::string Function_F_Id, Function_G_Id;
	std::string RestrictionOperator_Id;
	std::string InterpolationOperator_Id;
	std::string s;
	size_t N;
	size_t coarest_M;
	size_t nu1;
	size_t nu2;
	size_t itr_num;
	std::string Pro_Id;
	bool is_end = false;
	MultigridSolver<1>* MS_1;
	MultigridSolver<2>* MS_2;
	VectorFunction<1>* F_1d;
	VectorFunction<2>* F_2d;
	VectorFunction<1>* G_1d;
	VectorFunction<2>* G_2d;
	double max_norm;
	double E1, E2;
	double epsilon;
	std::string d = "../results";
	mkdir(d.c_str(), 0777);

	std::ifstream data_in("Inputfile");
	data_in >> c;
	while(!is_end) {
		if(c == '$') 
			std::getline(data_in, s, '$');
		data_in >> Pro_Id >> DIM >> N;
		data_in >> MultigridSolver_Id >> RestrictionOperator_Id >> InterpolationOperator_Id;
		data_in >> Function_F_Id >> Function_G_Id;
		data_in >> coarest_M >> nu1 >> nu2;
		std::string file_name = "../results/" + Function_F_Id + "_in_Problem" + Pro_Id + "_" + MultigridSolver_Id;
		file_name += "_" + RestrictionOperator_Id + "_" + InterpolationOperator_Id;
		size_t grid_num = 1 << N;
		file_name += "_" + std::to_string(grid_num) + ".txt";
		
		std::ofstream data_out(file_name);
		if(DIM == 2) {
			MS_2 = MultigridSolverFactory<2>::Instance()->CreateMethod(MultigridSolver_Id, RestrictionOperator_Id,
																		InterpolationOperator_Id);
			F_2d = VectorFunctionFactory<2>::Instance()->CreateFunction(Function_F_Id);
			G_2d = VectorFunctionFactory<2>::Instance()->CreateFunction(Function_G_Id);
		}
		else {
			MS_1 = MultigridSolverFactory<1>::Instance()->CreateMethod(MultigridSolver_Id, RestrictionOperator_Id,
																		InterpolationOperator_Id);
			F_1d = VectorFunctionFactory<1>::Instance()->CreateFunction(Function_F_Id);
			G_1d = VectorFunctionFactory<1>::Instance()->CreateFunction(Function_G_Id);
		}
		size_t l = 1;
		for(size_t i = 0; i != DIM; i++) {
			l *= (1 << N) + 1;
		}
		std::vector<double> v(l, 0);
		size_t num_step = 1;
		size_t cnt = 0;
		if(Pro_Id == "III")
			epsilon = 2.2*0.0000000000000001;
		else 
			epsilon = 0.00000001;
		if(Function_F_Id == "f1") 
			E2 = exp(sin(1));
		else if(Function_F_Id == "f2")
			E2 = 0.25;
		else if(Function_F_Id == "f3")
			E2 = exp(sin(1)*sin(1));
		else
			E2 = 0.25*0.25;
		data_out << std::fixed << std::setprecision(18);
		data_out << "grids" << '\t' << "num_step" << '\t' << "absE" << "\t\t" << "Reduction ratio" << "\t\t" << "Order p" << std::endl;
		data_out << (1<<N) <<'\t'<< 0 << '\t' << E2 << '\t' << std::endl;
		std::cout << Function_F_Id << "_in_Problem" << Pro_Id << " with " << MultigridSolver_Id << " with " << RestrictionOperator_Id;
		std::cout << " and " <<  InterpolationOperator_Id << " on grid with " << grid_num << "knots." << std::endl;
		std::cout << "step:" << cnt << '\t' << "error" << E2 << std::endl;
		size_t NUM = 16;
		if(Function_F_Id == "f3" && MultigridSolver_Id == "VC") 
			NUM = 6;
		if(Function_F_Id == "f4" && MultigridSolver_Id == "VC")
			NUM = 6;
		while(num_step < NUM && !( E2 <= epsilon && num_step >= 3)) {
			v.clear();
			v.resize(l);
			E1 = E2;
			if(DIM == 1) {
				(*MS_1).Solve(*F_1d, *G_1d, v, N, coarest_M, num_step, nu1, nu2, epsilon);
				E2 = (*MS_1).get_E(*G_1d);
			}
			else if(DIM == 2) {
				(*MS_2).Solve(*F_2d, *G_2d, v, N, coarest_M, num_step, nu1, nu2, epsilon);
				E2 = (*MS_2).get_E(*G_2d);
			}

			cnt++;
			std::cout << "step:" << cnt << '\t' << "error" << E2 << std::endl;
			data_out << (1<<N) << "&\t" << cnt << "&\t" << E2 << "&\t" << E2/E1;
			if(num_step > 0) {
				data_out << "&\t" << log2(E2/E1)/log2(1.0/(1<<N));
			}
			data_out << "\\\\" << std::endl;
			if(MultigridSolver_Id == "FMG" || MultigridSolver_Id == "FMG_1d") {
				num_step += 3;
			}
			num_step++;
		}
		data_out.close();
	//	std::getline(data_in, s, '%');
		data_in >> c;
		if(c == '#')
			is_end = true;
	}

	return 0;	
}
